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Abstract 



A great many observables seen in intermediate energy heavy ion collisions can be 
explained on the basis of statistical equilibrium. Calculations based on statistical 
equilibrium can be implemented in microcanonical ensemble (energy and number 
of particles in the system are kept fixed), canonical ensemble (temperature and 
number of particles are kept fixed) or grand canonical ensemble (fixed temperature 
and a variable number of particles but with an assigned average). This paper deals 
with calculations with canonical ensembles. A recursive relation developed recently 
allows calculations with arbitrary precision for many nuclear problems. Calculations 
are done to study the nature of phase transition in intermediate energy heavy ion 
collision, to study the caloric curves for nuclei and to explore the possibility of 
negative specific heat because of the finiteness of nuclear systems. The model can 
also be used for detailed calculations of other observables not connected with phase 
transitions, such as populations of selected isotopes in a heavy ion collision. 

The model also serves a pedagogical purpose. For the problems at hand, both 
the canonical and grand canonical solutions are obtainable with arbitrary accuracy 
hence we can compare the values of observables obtained from the canonical calcula- 
tions with those from the grand canonical. Sometimes, very interesting discrepancies 
are found. 

To illustrate the predictive power of the model, calculated observables are com- 
pared with data from the central collisions of Sn isotopes. 

Key words: heavy ion intermediate energy composites multiplicity 
PACS: 25.70.-z 25.75.Ld 25.10.Lx 
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1 Introduction 



Consider a central collision of two heavy ions. Nucleons from one nucleus will 
collide with nucleons from another nucleus. After a few collisions a given nu- 
cleon may lose the identity of its source. The system then more resembles a 
hot fluid of nucleons in an overall volume. Depending upon the original beam 
energy, this system may undergo an initial compression and then begins to de- 
compress. During this time the nucleons will interact with each other, at least 
between the nearest neighbours. As the density of the system decreases, higher 
density regions will develop into composites. As this collection of nucleons be- 
gin to move outwards, rearrangements, mass transfers, nucleon coalescence 
and most physics will continue to happen until the density decreases so much 
that the mean free paths for such processes become larger than the dimension 
of the system. Subsequently the objects follow the long range coulomb trajec- 
tories. Our objective is to have a soluble model which describes the physics of 
the situation at this freezeout density when one averages many nucleus-nucleus 
collisions. 

Although we chose central collisions to describe this scenario a similar situation 
will arise even for semi-central or semi-peripheral collisions. In such cases one 
may have projectile like fragment (and target like fragment and participants, 
region of violent collisions). For example, a projectile fragment may be excited 
which resembles a system of hot particles whose centre of mass velocity is close 
to that of the projectile [1]. 

The central assumption of the present article (and many others) is that equi- 
librium statistical mechanics can be used to describe the hot fluid of nucle- 
ons. Even the most well prepared experimental measurement of an energetic 
nucleus — nucleus collision represents an average of a very large number of 
initial states. In addition to this large number of different initial states, a 
large number of nucleon — nucleon collisions occur within each nucleus — 
nucleus collision. Together, this means that for many experimental observ- 
ables almost all the relevant phase-space can be opened up and described 
by the microcanonical ensemble in which the probability of reaching a chan- 
nel y is Q(y) / J2 y Q{y)- Here is the phase-space volume in the channel 
y. In the canonical ensemble the corresponding expression [2] is writen as: 
exp(—f(y)/T)/J2 ex P(—f(y)/T). Here f(y) is the free energy in the channel 
y. Since f(y) = —TlnQ{y) where Q(y) is the canonical partition function 
in the channel y, an equivalent expression is Q(y)/ J2Q(y)- A more detailed 
discussion of statistical equilibrium using reaction rates is given in Appendix 
A. 

The obvious experimental observables in heavy ion collisions are the number 
of nucleons and composites and their velocity distributions that result af- 
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ter the collision. The calculation of these in equilibrium statistical mechanics 
for Bevalac physics is more than twenty five years old [3,4,5]. At that time, 
the grand canonical ensemble was used to describe data from the Bevalac 
which normally used beam energies higher than 250 MeV/nucleon. However, 
at these energies most of the subtle and interesting features of equilibrium 
statistical mechanics as it pertains to heavy ion collision have disappeared. 
As the cross-sections of composites fall rapidly with A, the mass number, the 
most interesting results were productions of new particles such as pions and 
kaons, which can be included in the statistical model. Some discussion of this 
production is also given in Appendices A and C. Even so, Bevalac experiments 
brought out beautiful features of dynamics and established narrow limits on 
compressibility of nuclear matter and the momentum dependence of the real 
part of the optical potential. 

The applications of equilibrium statistical mechanics for intermediate energy 
heavy ion collisions started in the eighties. At these energies, the efforts 
switched to microcanonical ensembles although the concept of temperature 
was sometimes used [6,7,8]. One model called the Copenhagen SMM (statisti- 
cal multifragmentation model) is frequently used [8]. Another popular model 
is the Berlin model [7]. The use of the canonical ensemble, the main topic of 
this paper, is more recent [9]. It is as easy to implement as the grand canonical 
(and more accurate since fluctuations in the number of particles are elimi- 
nated: these sometimes cause large errors in computations of observables) . It 
is orders of magnitudes simpler than the microcanonical ensemble although 
in the latter more fine tuning can be done. These fine tunings do not appear 
important for most observables. 

What are the important issues we want to learn about in intermediate en- 
ergy heavy ion collisions? For many, it is to extract from data signals of a 
liquid-gas phase transition in nuclear matter. Nuclear matter is a hypotheti- 
cal large chunk of matter with N = Z where the coulomb interaction has been 
switched off. The p—V diagram for nuclear matter with reasonable forces looks 
like a Van der Waals equation of state [10]. One would then expect to see a 
liquid-gas phase transition if the experimental conditions are optimal. Such 
optimal conditions are discussed by Curtin, Toki and Scott [11] and Bertsch 
and Siemens [12]. For Bevalac energies the evolution of the temperature would 
go above the phase transition temperature but accelerators at the National 
Superconducting Cyclotron Laboratory (NSCL), the Texas A&M cyclotron, 
the Grand Accelerateur National D'ions Lourds (GANIL) and at Gesellschaft 
fur Schwerionenforschung mbH (GSI) can reach the liquid-gas phase transition 
region and offer the best possibility for experimental study. Further details of 
theoretical considerations which prompt an experimental investigation of the 
liquid-gas phase transition can be found in [13]. 

Unfortunately, this investigation of liquid-gas phase transition in intermediate 
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energy heavy ion collisions is fraught with many difficulties. Phase transitions 
occur in very large systems. In nuclear collisions, we are limited to three to 
four hundred nucleons (sometimes much less). For finite systems, signals of 
phase transition get diluted and distinctions between first order and second 
order transitions get blurred. The coulomb interaction, which prevents large 
nuclei from forming also interferes with the signals. It is thus necessary to 
use theories to clarify the situation. If one has a theory which fits many data, 
not necessarily related to phase transitions, but in addition predicts a phase 
transition one has some hope for the model to be valid. In this paper we will 
discuss phase transitions but in addition, data which will be compared to the 
thermodynamic model predictions. 



2 The Basic Formulae 



This section sets up the basic formulae of the model [9,14]. 

If there are A identical particles of only one kind in an enclosure at temperature 
T, the partition function of the system can be written as 



Qa = ^) A (1) 

Here oj is the partition function of one particle. For a spinless particle without 
any internal structure to = j^(2irmT) 3 / 2 ; m is the mass of the particle; V is 
the available volume within which each particle moves; A\ corrects for Gibb's 
paradox. If there are many species, the generalisation is 



Qa = T,U { ^t- (2) 

Here uii is the partition function of a composite which has i nucleons. For a 
dimer % = 2, for a trimer i = 3 etc. Eq. (2) is no longer trivial to calculate. The 
trouble is with the sum in the right hand side of eq. (2). The sum is restrictive. 
We need to consider only those partitions of the number A which satisfy 
A = X^rij. The number of partitions which satisfies the sum is enormous. 
We can call a given allowed partition to be a channel. The probablity of the 
occurrence of a given channel P{n) = P(ni, ri2, n^....) is 



m-^ii^ 1 (3) 

(4 a nil 
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The average number of composites of % nucleons is easily seen from the above 
equation to be 



(rii) = LOi 




(4) 



Since X^n^ = A, one readily arrives at a recursion relation [15] 



1 A 

k=l 



(5) 



For one kind of particle, Qa above is easily evaluated on a computer for A 
as large as 3000 in matter of seconds. It is this recursion relation that makes 
the computation so easy in the model. Of course, once one has the partition 
function all relevant thermodynamic quantities can be computed. 

We now need an expression for uo k which can mimic the nuclear physics situ- 
ation. We take 



where the first part arises from the centre of mass motion of the composite 
which has k nucleons and q k is the internal partition function. For k — 1, 
q k = 1 and for k > 2 it is taken from the Fermi-gas model. For each composite 
consisting of k nucleons, we approximate the intrinsic free energy at freeze-out 
by E - TS = -W k + a(T)k 2 / 3 + kT 2 /e - T x 2kT/e where e is a constant. 
This gives 



Here, as in [8], W / o=16 MeV is the volume energy term, cr(T) is a temperature 
dependent surface tension term. The value of eo is taken to be 16 MeV. The 
explicit expression for a{T) used here, as in [8], is 



with o"o =18 MeV and T c = 18 MeV. In the nuclear case one might be tempted 
to interpret V of eq.(6) as simply the freeze-out volume but it is clearly less 
than that; V is the volume available to the particles for the centre of mass 
motion. Assume that the only interaction between clusters is that they can 
not overlap one another. Then in the Van der Waals spirit we take V = 
Vf r eeze—V ex where V ex is taken here to be constant and equal to Vq = A/ p . The 
assumption that the interaction between different composites is only reflected 



(6) 



q k = exp[{W k - a{T)k 2 / 3 + T 2 k/e )/T] 



(7) 



a(T) = a [(T 2 - T 2 )/(T 2 + T 2 )] 5 / 4 
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through an excluded volume and that this excluded volume is independent 
of multiplicity is an idealisation which will fail for a non-dilute system. We 
therefore restrict the model, somewhat arbitrarily to volumes Vf reeze > 2V . 
There are experimental signatures that Vf reeze is indeed greater than 2V [13] 
so this is not a debilitating feature of the model. In all our considerations we 
restrict p/p to less than 0.5. 

Among quantities of interest is the inclusive cross-section given by eq.(4). Ac- 
tually this is a simplification. The occupation given by eq.(4) is the occupation 
of the composite with i nucleons at temperature T. Both the ground state and 
the excited states contribute to (rii). Some of the excited states will be parti- 
cle unstable and will decay into lower mass composites before they reach the 
detector. On the other hand, some higher mass composites, will, by the same 
argument, decay into the composite i. In later sections, where we compare 
populations with data, this aspect will be taken care of. The expression for E 
at a given temperature T is simple (this is needed for a caloric curve which is 
measured in experiments). The energy carried by one composite is given by 

E k = T 2 dlnuo k /dT = \T + k(-W + T 2 /e ) + a(T)k 2 / 3 - T[da{T) /dT}k 2 ' 3 . 

Of these the first term comes from the centre of mass motion and the rest from 
q%. The term T[da{T) / dT]k 2 / 3 comes from the temperature dependence of the 
surface tension. It has a small effect. The energy of the whole system is given by 
E = T 2 ^^§f. Using eqs.(2) and (4) we arrive at a very transparent formula 
:E = J2( n k)Ek- The pressure is given by p = TdlnQA/dV. If for purposes of 
illustration, we neglect the long range Coulomb interactions and use eqs.(2) 
and (4) we get p = Ty J2( n i)- This is just the law of partial pressures. 

For analysing phase transitions in the model, it is very useful to calculate the 
average value of the largest cluster in the ensemble. Eq. (2) shows that the 
size of the largest cluster varies. In that ensemble there is a term For 
this the largest cluster is the monomer. For example in eq.(2) we also have 

n (A/2-n/2) 

a term (A/2-n/2)! • Here the largest cluster is the dimer. Consider building 

Qa with u>i, UJ2 uJk, 0, 0, 0, 0... In this ensemble the largest cluster will span 

from monomer upto a composite with k nucleons. Let us label this partition 
function Qa(^i, ^2, ■■■■^k, 0, 0, 0..). Let us also build a Qa where the largest 
non-zero u> is u>k-i- The partition function is Qa(wi, <^2, ••k'fc-i, 0, 0, 0, 0). In 
this ensemble all the previous channels are included except where the largest 
cluster had k nucleons. If we define 



AQaO) = QA(vi,V2-Vk, 0, 0, ..) - Qa(vi, u 2 ...Uk-i, 0, 0....), 
then the probability of the largest cluster having k nucleons is 
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Pr(k) 



&Q A (k) 
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Qa(uji,uj 2 uj a ) 



If we now label the average value of the largest cluster as (k max ), then (k max ) = 
Y,k x Pr{k). A more useful quantity is i^pl. The limits of this are and 
1. 

Another interesting quantity which has been the subject of an enormous 
amount of interpretation [17] is the multiplicity distribution of a species or 
a group of species. In most models this requires a very elaborate Monte-Carlo 
calculation. In the canonical ensemble there is an elegant equation. 

1 ui n 

P n{k) = -pz jQA-nk(vi, W 2 — Wfc-1, ^fc = 0, U k+1 ..U A ) (9) 

Here P n (k) is the probability of obtaining the composite k n times. 

The strength of the canonical model as described here lies in the fact that 
all calculations above avoid Monte-Carlo sampling. In many other models, a 
Monte-Carlo sampling over the channels is required. Since the number of chan- 
nels is enormous, this requires great ingenuity but also much more computer 
time. 

The model of one kind of particles in which composites have a volume energy, 
a surface energy and excited states is already very useful for investigations of 
phase transition, caloric curves etc. and we will pursue this in latter sections 
a great deal. Let us nonetheless introduce here the model with two kinds of 
particles (so that one can compare with actual nuclear cases) [13,14,16]. Now 
a composite is labelled by two indices u — > Uij. The partition function for a 
system with Z protons and N neutrons is given by 



There are two constraints: Z = J^i x n i,j an d N = J2j x n i,j- These lead to 
two recursion relations any one of which can be used. For example 




(10) 



Qz,N - 



ivijQ Z -i,N-j 



(11) 



where 



V 



(2-nmTfl\i + jf' 2 x q,, 



(12) 
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Here qij is the internal partition function. These could be taken from experi- 
mental binding energies, excited states and some model for the continuum or 
from the liquid drop model in combination with other models. The versatility 
of the method lies in being able to accommodate any choice for q i:j . A choice 
of qij from a combination of the liquid drop model for binding energies and 
the Fermi-gas model for excited states that has been used is 



q id = exp -[W Q a - aa 2/s - - s± ^- + T 2 a/e ] (13) 

where a = i + j, W = 15.8 MeV, a = 18.0 MeV, k = 0.72 MeV, s = 23.5 MeV 
and 6q = 16 MeV. One can recognise in the parametrisation above, the volume 
term, the surface tension term, the coulomb energy term, the symmetry energy 
term and contributions from excited states. 



The coulomb interaction is long range. Some effects of the coulomb interaction 
between different composites can be included in an approximation called the 
Wigner-Seitz approximation. We assume, as usual, that the break up into 
different composites occurs at a radius R c , which is greater than the normal 
radius Rq. Considering this as a process in which a uniform dilute charge 
distribution within radius R c collapses successively into denser blobs of proper 
radius R it j, we write the coulomb energy as [8] 



Ec= i^r + ^5R- (1 -^ w (14) 

It is seen that the expression is correct in two extreme limits: very large freeze- 
out volume (R c — > oo) or if the freeze out volume is the normal nuclear volume 
so that one has just one nucleus with the proper radius. 

For the thermodynamic model that we have been pursuing, the constant term 
f^ 5 - is of no significance since the freeze-out volume is assumed to be con- 
stant. In a mean-field sense then one would just replace the coulomb term in 
eq.(13)by^(1.0-(p/p ) 1/3 ). 

Before we leave this section, we mention that the mass parametrisation im- 
plied by eq.(13) can be vastly improved with only slight complications. We 
will later present results with the improved formula [16]. A pedagogical issue: 
although we have derived results here based on eq.(l) which takes care of 
(anti)symmetrisation only approximately it is shown in [18] that the specific 
structure of eqs.(5) and (11) occur more generally when (anti)symmetrisation 
is included properly. Part of this argument is presented in appendix B which 
also demonstrates that results based on this section are quite accurate. 
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3 General features of yields of composites 

We pursue here the model of one kind of particles. For 200 particles at a 
constant freeze-out volume =3.7Vo we have plotted in fig.l (n k ) (in the figure 
we call this Y (a)=yield of composite of mass a) at three temperatures. At the 
lowest temperature shown the curve has a U shape. The yields Y (a) first begin 
to fall, then reach a minimum and then the yields for heavier masses increase 
finally cutting off at 200. In the literature the heavy fragments are called 
the liquid phase. The light fragments are gas particles. As the temperature 
increases the maximum at the higher a decreases in height finally disappearing 
at ~6.35 MeV. At higher temperature Y(a) falls monotonically. The surface 
tension plays a crucial role in this evolution. At any temperature the lowest 
value of the free energy E — TS will be obtained. It costs in the energy term E 
to break up a system. A nucleus of A nucleons has less surface than the total 
surface of two nuclei each of A/2 nucleons (the volume energy term has no 
preference between the two alternatives). Therefore at low temperature one 
will see a large chunk. The -TS term favours break up into small objects. The 
competition between these two effects leads to the general features of Y(a) as 
a function of temperature. As we will see in the next section, the temperature 
at which the maximum of the yield at the high side of a disappears is the 
phase transition temperature. 

Similar features are seen also in other models of multifragmentation as applied 
to nuclear physics. The earliest such model was the percolation model [19,20]. 
The model has a parameter p which gives the probability of two nearest neigh- 
bour sites joining together as in a composite. Beyond a certain value of p, a 
percolating cluster is formed which goes from edge to edge of the system. This 
corresponds to the large cluster which forms at the lower temperature in fig.l. 
The lattice gas model [21] has similarity with the percolation model but has 
a Hamiltonian, includes percolation model as a subset [22] and also includes 
the formation of a percolating cluster 



4 Phase transition in the model 

4-1 Signatures from thermodynamic variables 

We now begin the discussion of a phase transition in the model. The free energy 
of a system of particles is given by F = — T\u.Qa and IiiQa is directly calculable 
with eq.(5). For a system of 200 and 2000 particles, the free energy per particle 
is shown in the top panel of fig. 2, as a function of temperature for fixed freeze- 
out density 0.27p - An approximate break in the first derivative of F/A is 
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seen to develop at ^6.35 MeV for 200 particles and at «7.15 MeV for 2000 
particles. We believe the break would be rigorous if we could go to an infinite 
system. A break in the first derivative implies a first order phase transition and 
a discontinuous change in the value of entropy per particle. This would imply 
that the specific heat at constant volume per particle c v = (^ffr^)y would go 
through a peak (for an infinite system this peak would go to oo). We show this 
in the midddle panel of fig. 2 for systems of 200 and 2000 particles, where we 
find the width of the peak decreases and the height of the peak increases as 
the particle number increases. As expected, the temperature where the specific 
heat maximises also coincides with the temperature at which the maximum 
in the high side of a (fig.l) just disappears. 



10' 



10' 



10 3 



10 s 



10" 




T = 5.80 MeV 
T = 6.35 MeV 
T = 7.30 MeV 



Fig.l Yield Y(a) against a at three different temperatures. The dissociating 
system has 200 particles and the freeze-out density is 0.27po- 



Another very interesting quantity is the quantity (k max )/A (i.e. the size of 
the largest cluster) as the temperature varies. This can be calculated using 
eq.(8).We define as the temperature where the break in the derivative of 
the free energy occurs (this is the first order phase transition temperature). 
Calculating the size of the largest cluster at different temperatures, we find 
that (k max )/A approaches 1 as T < T b and approaches a small number as 
T > Tjj. The change is smooth for low mass nuclei (bottom panel of fig. 2) but 
becomes more sudden for larger systems. For large systems there is a large 
blob (i.e., liquid) below T b which disappears as soon as T crosses T b . This 
we think is a very engaging example of boiling emerging from a theoretical 
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calculation. 
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Fig. 2. The free-energy per particle, the specific heat at constant volume, Cy j A 
and the size of the largest cluster as a function of temperature, systems of 200 
and 2000 particles. 

To summarise, the thermodynamic model predicts unequivocally a first order 
phase transition at intermediate energy. In the realm of density pj po < 0.5 for 
which we believe the model to be reasonable there is no critical point. Bugaev 
et al. [23] have taken the model beyond this range of density and find that the 
critical density is p/po — 1 and the temperature is 18 MeV when the surface 
tension a(T) goes to zero. We end this section by noting that microcanonical 
calculations using statistical equilibrium were also suggestive of a first order 
phase transition occurring at intermediate energy [24,6]. 



4-2 Power-law and scaling behaviour of composite yields 

A rather large part of literature in heavy ion reaction postulates that in mul- 
tifragmentation at intermediate energy one is near the critical point of nuclear 
matter. One then proceeds to determine from the data the critical temperature 
and various critical exponents. The working formula, obtained from models of 
critical phenomena (to see how the formula arises, refer to [25,26]) is 
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(n a ) = a- T f(a°(T-T c )) 



(15) 



Here r is called the Fisher exponent [27], a is the mass number of the com- 
posite, a is a critical exponent and T c the critical temperature ; / is as yet an 
unspecified function but instead of being a general function of a and T it is a 
function only of the combination a a (T — T c ). This is called scaling. At T = T c 
the yield (n a ) = a~ T /(0) is a pure power law but away from T c it will deviate 
from a power law. 

In intermediate energy collisions, even if we proceed under the assumption 
that one is observing critical phenomena we can not expect near perfect fit to 
eq.(15) whose validity depends upon the dissociating system being very large. 
Also the range of a is to be chosen judiciously. It can not be very small (since 
eq.(15) applies to "large" a's [25]). But a also should be truncated on the high 
side (significantly smaller than the size of the dissociating system). 

With these provisos we can at best expect a moderately good fit. Extracting 
r, a and T c from a given set of (n a ) (either from experiment or models) when 
only an approximate fit is expected is non-trivial and not unique. We skip 
the details here which are given in [28,29,30]. A more sophisticated method of 
extraction of the relevant parameters can be found in [31]. The same technique 
is used in [32]. 

The EOS collaboration [33] obtained data from break up of 1.0 GeV per 
nucleon gold nuclei on a carbon target. Depending upon the impact parameter, 
the excitation energy (or the temperature) of the projectile like fragment which 
breaks up into many composites will vary. In [28,29] it is argued that T in 
eq.(15) varies linearly with the charged multiplicity m and the scaling function 
of eq.(15) is changed from f(a a (T—T c )) to /(a°"( m ~ m-c )). Here m c is the critical 
multiplicity. Having determined from the data r, a and m c (as mentioned 
before we are skipping the details of how the extraction is done but this 
can be found in [28,29]) one then verifies if the scaling law works: that is, 
we check if for all a's, {n a )a T will fall on the same curve when plotted as a 
function of a u [m — m c )/m c . How well this works can be seen, for example, in 
fig. (18) of [30]. The deviations from the hypothetical "universal" curve are by 
no means negligible but can we assume that the scatter of points is entirely 
finite particle number effect and conclude that we have indeed seen evidence 
of critical phenomena? 

To resolve this, we play a theoretical game. We take the thermodynamic model 
(which we know has only a first order phase transition), pick a system with 
particle number A, generate (n a ) for different temperatures T and from these 
data extract best possible values of r, a and T c . Having obtained these we 
examine how well the scaling law applies. This is shown in fig. 3. The figure 
is taken from [32] where other similar examples are displayed. For the model, 
deviations from one "universal" curve are smaller than what the EOS collab- 
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oration data gave. We might conclude we have extracted the model critical 
temperature and the critical exponents. These would be wrong conclusions of 
course because the model has only a first order phase transition. In fact the 
value of T c one extracts this way is quite close to T b , the first order phase 
transition temperature. 




a (T-T c ) 



Fig. 3. The scaling behaviour in the mass range 10 < a < 40 in the thermody- 
namic model for different systems at different freeze-out densities. 



In [30] the Copenhagen SMM is used to show that approximate scaling is ob- 
tained. The hope would then be that the theory also demonstrated criticality. 
The SMM is, in spirit, very close to the thermodynamic model, thus we doubt 
that the very approximate collapse of (n a )a T on one curve is any indication of 
criticality. It is impossible to disentangle what errors arise because the wrong 
formula is applied and what errors arise because of the finiteness of the system 
and many other factors such as the coulomb force, pre-equilibrium emission 
etc. Experimental data would have a very hard time of choosing between a 
first and a second order transition. 
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From 8 GeV/c 7r~ on Au data, the ISiS [34] collaboration obtained the caloric 
curve [35]. The specific heat was obtained by differentiating with respect to 
T. Experiment shows that the peak of the specific heat coincides well with the 
position where the x 2 f° r J2(( n a) — Ca~ T ) 2 minimises. Here both C and r are 
taken as parameters to be fixed by minimisation. The canonical model gives 
similar results. Further details of experiment and theory can be found in [36] 
where effects of the coulomb interaction on the position of the maximum of 
the specific heat is discussed in detail. 

We turn now briefly to another phenomenological model which was invoked 
twenty years ago [37] but was revived recently [38]. This is yet another example 
where evidence for criticality can be drawn too hastily. Consider the formation 
of a droplet containing a particles in the liquid phase surrounded by b particles 
in the gas phase. At constant temperature and pressure the Gibbs free energy 
is the relevant factor. Then 

G wit hdrop = Via + figb + 4irR 2 a + Tt lna 

and G no drop = H g (a + b) 

The probability of formation of droplet containing a particles is proportional 
to exp(— AG/T) so that the yield of droplets of size a is 



{n n 



Ca~ T exp[(/i s - ^)a/T + c 2 a 2/3 /T] (16) 



Here both \ii and /i g are functions of T. At coexistence and also at critical 
temperature, they become equal to one another. Also c 2 is a function of tem- 
perature and at T c , the coefficient c 2 goes to zero. Since above T c , there is 
no distinction between the liquid and the gas phase, one can not speak of 
droplets. Thus the theory only applies to T < T c . As such, the formulation is 
more limited than that of eq.(15) which applies to both sides of T c . We now 
generate values of (n a ) from the thermodynamic model for different temper- 
atures and try to fit these "data" using eq.(16). The following fit was tried. 
We set t—2. Let a = (fj, g — fJ>i)/T, 7 = c 2 /T. We fit the calculated (n a ) to 
Ca~ 2 exp(cra + 7a 2//3 ) at different temperatures where a, 7 values at each tem- 
perature are varied for best fit. The values of a, 7 as functions of temperature 
are shown in fig. 4 where we also show that the parametrisation fits the values 
of (n a ) very accurately. The values of a and 7 both go to zero near temper- 
ature T=6.5 MeV suggesting the critical temperature is 6.5 MeV. Of course 
this conclusion would again be wrong since the model which gave these (n a ) 's 
has only first order phase transition. 

One problem is that whenever a fit, whether through eq.(15) or through 
eq.(16), is done, the fit is attempted for a narrow range, a=6 to 40. In this 
limited range moderate to excellent fits are obtainable for different looking 
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parametrisations. It is shown in [32] that if the range of a could be extended 
to beyond 100, different parametrisations would diverge. Unfortunately, the 
range of a has to be limited. For example higher values of a would have contam- 
ination from fission processes which is something we do not wish to include. 
If we are stuck to a limited range of a's we will also be limited by ambiguity. 
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Fig. 4. The parameters of the droplet model a and 7 as a function of temper- 
ature for a system of 240 particles at freeze-out volume 4Vo- The right panels 
show the fit of the model to the yields obtained in the thermodynamic model 



The emphasis towards unravelling critical phenomena from data on intermedi- 
ate energy heavy ion collisions is at least partly due to history. The observation 
by the Purdue group [39] that the yields of the fragments produced in p + Xe 
and p + Kr obeyed a power law (n a ) oc a" T led to a conjecture that the frag- 
menting target was near the critical point of liquid-gas phase transition. The 
origin of this conjecture is the Fisher model [27] which predicts that at the 
critical point the yields of the droplets will be given by a power law. Also the 
first microscopic model that was used [19,20] to compute yields of fragments 
was the percolation model which has only a continuous phase transition and a 
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power law at criticality. The power law is no longer taken as a "proof" of crit- 
icality. There are many systems which exhibit a power law: mass distribution 
of asteroids in the solar system, debris from the crushing of basalt pellets [40] 
and the fragmentation of frozen potatoes [41]. In fact the lattice gas model 
which has been used a great deal for multifragmentation in nuclei gives a 
power law at the critical point, at the co-existence curve (this is a first order 
phase transition provided the freeze-out density is less than half the normal 
density) and also along a line in the T — p plane away from the coexistence 
curve [42,43,44]. 

We finish this section stating that the lattice gas model which has a Hamilto- 
nian and can be and has been used to fit many data (not in any obvious way 
connected with phase transition) also predicts a first order phase transition at 
intermediate energy [42,43]. 

4-3 Comparision with mean- field theory 

Here we concentrate on the thermodynamic model but as applied to nuclei 
with neutrons and protons. The operative equations are (10) to (13) but we 
will switch off the coulomb term (k of eq.13 will be set to zero). The objective 
is to compare with finite temperature Hartee-Fock results for nuclear matter. 
For nuclear matter the coulomb interaction has to be switched off and one 
retains only the nuclear part of the interaction. 

Phase transitions are often considered in the mean-field model. Examples for 
the present discussion are [10,12,13]. Invariably a grand canonical ensemble 
is used characterised by a neutron chemical potential p^(T,p) and a proton 
chemical potential pp(T, p). The use of the grand canonical model would imply 
that the results are valid for very large systems although in nuclear physics 
we often use the grand canonical ensemble for not so large systems as well. 

Muller and Serot [45,46] used the mean- field model to investigate phase tran- 
sition in nuclear matter. Normally nuclear matter means a very large system 
with N = Z with the coulomb force switched off. For this section we will use 
the term nuclear matter for very large systems but N can be different from Z. 
The coulomb is switched off as usual. Define proton fraction y = Zj (N + Z). 
Symmetric nuclear matter has y=0.5 and would have a first order phase tran- 
sition below the critical point. But for y deviating significantly from 0.5, these 
authors demonstrate with a more general Maxwell like construction that the 
first order phase transition would turn into second order. Further the phase 
transition would take place neither at constant volume nor at constant pres- 
sure but would have a more general path to traverse. 

The general characteristics of mean-field theories is that one is constrained to 
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have one density. Having the same density everywhere is a big price to pay. For 
example, this would not permit a liquid phase at one place and a gas phase at 
another. The limitation of one density only shows up as mechanical instability, 
i.e., in parts of the equation of state diagram (p — p isothermals) dp/ dp turns 
out to be negative. This is unacceptable for infinite matter and then one has 
to, by hand, correct this using a Maxwell construction. The thermdynamic 
model is very different. Here, for example, p/po = 0.3 does not mean that at 
the freeze-out volume, matter is uniformly stretched. Rather matter breaks up 
into different blobs all with the same density p but there are empty spaces 
between blobs. If there is a large blob, we identify it as liquid, nucleons and 
light composites in the adjoining spaces form the gas (In [47], it is shown 
that this last scenario has a lower free energy compared to uniform stretching 
as asumed in Hartree-Fock theory). For large matter, there is no need for 
(dp/dp) T to be negative. 

A similar thing happens with isospin fractionation. In mean-field theory, there 
is one value of y everywhere. Experimentally, it is verified that if the dissoci- 
ating system has a small y, then after break up, the largest blob has y > y<n ss 
whereas n p j {n p + n n ) < y diss where Tip, 77^ £1X6 free protons and free neutrons 
respectively. Here y^ss is the the y value of the dissociating system. One might 
say the liquid phase has a different y from that of the gas phase. Again mean- 
field theory would have a hard time accommodating this. It must have the 
same value of y everywhere and the fact that this is an unstable situation shows 
up in the following way. If we draw Pp(pn) as a function of y at constant tem- 
perature, the derivative (d/ip/dy) p can turn out to be negative (equivalently 
(dpN /dy) p can turn out to be positive). In the thermodynamic model, isospin 
fractionation happens naturally. In general, the model has, as final products, 
all allowed composites, a,b,c,d...., where the composite a has y a = % a j {ia+ja) 
where i a ,j a are the proton and neutron numbers of the composite a. The only 
law of conservation is Z = J2 a ^a x n a and N = J2 a Ja x n a - So a large chunk 
can exist with higher y than that of the whole system and populations of other 
species can adjust to obey overall conservation laws. Whatever partition low- 
ers the free energy will happen. Since we are using a canonical model, we do 
not need the chemical potentials \ip or p N but we can compute them anyway 
from the relation p = (dF/dn) v ,T- We know the values of Qz,n,Qz-i,n and 
Qz,n-i- Since F = —TlaQ, one has pp = — T (InQ Zt N-^Q z-i,n) and similarly 
for Pn • 

Calculations with the canonical model discussed in this article, do not show 
regions of negative (dp P /dy) P:T [48,49]. These also suggest that the phase 
transition in the canonical model remains first order for asymmetric matter. 
We show in fig. 5 results of cy calculation for different degrees of asymmetry. 
One sees that as the system gets bigger, the maximum in cy becomes narrower 
and higher, ensuring there will be a break in the first derivative of the free 
energy in the large matter limit. 
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Chemical instability for finite systems in Hartee-Fock theory has also been 
worked out. Contributions of both coulomb and surface terms can be included. 
For details see [51,52,53]. 
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Fig. 5. The CyjA as a function of temperature for systems of 200,500 and 1000 
particles with different proton fractions: (y = Z/A). 



5 Comparison of canonical and grand canonical 

As noted in the introduction, the grand canonical version of the model we are 
pursuing in this paper has been known and used for a long time. Now that 
we know how to treat an exact number of particles rather than an ensemble 
of particle numbers, it will be useful in a few cases to examine, given that our 
dissociating system has an exact number of particles, how the use of grand 
canonical ensemble affects the prediction of observables. For simplicity, we 
start with the model of one kind of particles and our dissociating system has 
200 particles. Thus we can have monomers, dimers, trimers,...upto a compos- 
ite of 200 particles. In the grand canonical ensemble, the average number of 
composites with k nucleons is given by 

(n k ) = exp(/3fj,k)u k = exp(Pfj,k)Vu k (17) 
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Here (3 is the inverse of temperature and uou is the same as defined in eq.6 
and \i is the chemical potential. We also use u = u/V where Cj only depends 
upon the composite and the temperature but not upon the volume of the 
dissociating system. The chemical potential is determined by solving 



p = ^kexp(k(3[i)u) k (18) 
k=i 

In this example k m =the number of particles in the largest cluster=200. Having 
determined /i we now find (n^) from {n k ) = exp(kj3fi)Vujk- 




Fig. 6. Comparison of yields obtained in the canonical and grand canonical 
models at different temperatures, for a system of 200 particles at freeze-out 
volume 3Vq. 



In figs. 6 and 7, we make a comparison of (nfc)'s from canonical and grand 
canonical ensembles. The value of V was set at 3Vo. Results are shown for 
temperatures below the phase transition temperature and above it. Fig. 6 
seems very reasonable. The overall features are similar. The differences get 
highlighted in fig. 7. At temperature 7.3 MeV, {n k ) GC and {rik) c are practically 
the same upto fc=40 but deviate wildly afterwards. Since most of the time we 
are not interested in the heavier products and k=40 is the limit of intermediate 
mass fragments one is investigating, the grand canonical ensemble does an 
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adequate job. We must be aware however, that, for heavy composites the grand 
canonical ensemble does a very poor job. The accuracy of the grand canonical 
ensemble at temperature 5.8 MeV (below the phase transition temperature) is 
absolutely awful for almost all composites. This is also the temperature range 
appropriate for most intermediate energy reactions. It is thus dangerous to 
use the grand canonical ensemble in intermediate energy heavy ion reactions. 
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Fig. 7 The ratio of yields obtained in the grandcanonical and the canonical 
model at different temperatures. 



If however, one is only interested in finding the ratio of populations of two 
adjacent composites, the grand canonical continues to be useful over a larger 
domain. This is shown in fig. 8. 




Fig. 8 Ratio of yields of adjacent composites in the two models. 
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The very different populations of composites below the phase transition tem- 
perature leads to drastically different caloric curves in the grand canonical 
ensemble and canonical ensemble. As noted in section 4 and shown in fig. 2, 
for a fixed density the specific heat per particle maximises at a certain temper- 
ature. Keeping density fixed, if we increase the number of particles the height 
of the maximum increases and the width decreases. In fig. 9 we show this again 
for 200 and 2000 particles, but now we have also indicated the specific heat cal- 
culated in the grand canonical ensemble. In both the models, the peak of the 
specific heat increases when we go from 200 to 2000 particles and the widths 
decrease but the results are much more dramatic in the canonical model. In 
particular, it is not obvious that the specific heat in the grand canonical en- 
semble will attain extraordinary heights and miniscule widths. In fact, it was 
suggested in the literature, engaging the grand canonical ensemble, that there 
is a discontinuity in the value of the specific heat at phase transition but no 
infinity [23]. 
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Fig. 9 Specific heat per particle at constant volume when the system has total 
number of particles 200 and 2000. Canonical and grand canonical values are 
shown. 

To understand at a more fundamental level the cause of the difference in values 
of specific heats in the two ensembles, we will analyse the case of 2000 particles 
in some detail. In the grand canonical model, even though we are using the 
average value of the particle number to be 2000, there are, in practice, systems 
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with varying particle numbers (in principle, from to oo). The part which has, 
for example, 1000 particles has density half of the prescribed density. The peak 
in the specific heat of this half density will occur at a different temperature 
than that which maximises the specific heat at density 2000/ V. Thus there 
is a smearing effect. This is always an inherent problem with using the grand 
canonical ensemble but most of the time the fluctuation from the average 
value is small enough that one can live with it. This would have meant, in our 
present example, the part which contains 1000 particles is so negligibly small 
that it does not matter. This however is not so in the present model below 
the phase transition temperature. 

In the present case, the grand canonical calculation starts out by obtaining p 
from eq.18 where /c m =2000; p was taken to be po/2.7. The average value of 
(rik) is then given by eq.(17) where V = 2000 x 2.7 /p . With this we have 
Efc=?°°° k{n k ) = 2000. The fluctuations in the model can be calculated easily. 
We have the general statistical relation 



1 d 2 lnQ gr ., 

p 2 Wp~ 



= (N 2 ) - (N) 2 (19) 



Here Q gr . C an is the grand canonical partition function. We can write two ex- 
pressions for Q gr.can- One is: 



fc=2000 

InQ gr.can = exp((3pk)u k (20) 
k=l 

This immediately leads to 



fc=2000 

(N 2 )-(N) 2 = £ k 2 (n k ) (21) 
k=i 

which is easily calculable. The other expression we can exploit in the present 
case is 



Qgr.can = eX P(^ K )Q K,k m (22) 

K=l 

where QK,k m is the canonical partition function of K nucleons but with the 
restriction that the largest cluster can not have more than k m {= 2000) nucle- 
ons. We can calculate these explicitly using methods of section 2. For practical 
reasons, K has to be cut off at the upper end. Here we used fT=10000 as the 
upper limit. Since the average number of particles is 2000, this appears to be 
a safe upper limit in eq. 22. The quantity p is known from solving eq. 18. 
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The fluctuations calculated with eqs. 21 and 22 are shown in fig. 10. One 
sees there is a temperature above which the fluctuations are small. At these 
temperatures, the grand canonical value of specific heat is indistinguishable 
from the canonical value. But as the temperature is lowered, fluctuations grow 
rapidly and the results begin to diverge. 
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Fig. 10 Fluctuations calculated using eqs. (21) and (22). The solid line corre- 
sponds to using eq. (22) with K cut off at 10,000 and the dotted line cor- 
responds to using eq. (21). T\ corresponds to to the temperature where the 
specific heat maximises in the canonical calculation and T<i to the temperature 
of highest specific heat in the grand canonical calculation. 

It is interesting to study fluctuations further. The probability of K parti- 
cles being in the grand canonical ensemble is oc e K P^+ ln QK _ We plot in fig. 11 
exp[(3n(K — A) + ItiQk — IuQa\. This takes the value 1 at K = A and in 
the normal picture of the grand canonical ensemble would drop off rapidly on 
either side of A. This does happen at temperatures higher than the boiling 
temperature. The case at T=7.7 MeV corresponds to a standard scenario. But 
the situation at temperature 7.3 MeV is drastically different. The probability 
does not maximise at K = A but at a lower value. It is also very spread out 
with a periodic structure. The periodicity is 2000 and is linked with the fact 
that in the case studied, the largest composite has 2000 nucleons and at low 
temperatures, this composite will play a significant role. 

More discussions on this case can be found in [50]. 
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Fig. 11. These graphs show the spread of particle numbers in the grand canon- 
ical ensemble when the average particle number is 2000. The spread is very 
narrow at temperature 7.7 MeV but becomes quite large at lower tempera- 
tures. 



6 Specific heat at constant pressure 



We have used Cy, the specific heat at constant volume a great deal in the 
previous sections. In canonical models Cy is always positive. Writing: (E) = 

^E^P^P and Cv = ^ E )/ dT ^ we § et °v = M(E~ (E)) 2 ) which 
is the expectation value of a positive definite operator. However, specific heat 
at constant pressure allows no such generalisations. Here we enter into a dis- 
cussion of the specific heat at constant pressure in the thermodynamic model. 
We should add that dissociation after two heavy ions collide is largely an 
uncontrollable situation and we do not know what is a better description: 
disassembly at constant volume, disassembly at constant pressure or a hybrid 
situation. 

Lately, interest in the topic has increased with the realisation that for finite 
systems, C p can sometimes be negative and such cases might arise in heavy 
ion collisions [54,55,56]. To study this possibility in our model, we find it 
convenient to look at the p—p diagram at constant temperatures (isothermals) . 
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This is shown in fig. 12. We see there are regions of mechanical instability where 
(fp) T < 0- We will show that the occurrence of negative C p happens in this 
region. 
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Fig. 12. EOS in the canonical model for a system of A=200 particles. The 
largest cluster also has iV=200. 

The most famous case of mechanical instability is the Van der Waals equation 
of state. In nuclear physics, if one uses the Hartree-Fock theory, then also 
large regions of mechanical instability appear. Examples of this can be seen 
in many published works: [10,9,48]. All these published works are for infinite 
systems (unlike the p — p diagram for fig. 12 which is drawn for 200 particles). 
Quantitative examination of the equation of state diagrams reveal that the 
regions of mechanical instability are far bigger in the case of Hartree-Fock 
as opposed to what we see in fig. 12. In fact, plotted on the same scale, the 
region of mechanical instability would be tiny (ref. fig. 1 of [48]) and one would 
have to plot it in an expanded scale (such as is done in fig. 12) to study it 
quantitatively. 

In the Van der Waals case or in the Hartree-Fock case for infinite nuclear 
matter one uses a Maxwell construction to replace the region of mechanical 
instability [2]. In the thermodynamic limit, regions of mechanical instability 
should disappear. In our case there is no prescription for Maxwell construction. 
Also since our system is very finite, we take the mechanical instability in fig. 12 
as real and follow the consequences for the specific heat. In the figure we have 



27 



drawn isothermals at three temperatures; T\ < T2 < T3. Here T2 is only slightly 
higher than T\. Instead of p let us use the variable V oc 1/p. The pressure is 
given by p = T{m/V) where m is the multiplicity. [We actually use m — 1 
but this is inconsequential for the discussion to follow.] For the simple case of 
monomers only, p is given by p = T{A/V) where A is the number of particles. 
This number does not change thus p keeps falling with V . In our case, m is 
significantly less than A. It is not a constant as V and/or T change. As can 
be readily guessed, m increases with T at constant V; m also increases with 
V at constant T. Negative compressibility is marked by (dm/dV) T > m/V. 



Let us consider the points c and d in fig. 12. Let c have multiplicity m, volume 
V and temperature T; for d the corresponding quantities are m + 5m, V + SV 
and T + 5T. Here 5V is negative, 5T is positive. Using 



p = T™ = (T + 5T) m + 6m (23) 

y V v J v + sv K J 



we arrive at 



5rn = 5V_5T 

m V T K J 



In the region (c, d), 5V is negative, ST is positive thus 5m is negative. If 
m goes down then so does the potential energy (creating more m creates 
more surface and hence increases energy). The change in kinetic energy is : 
|[(m + 5m)(T + 5T) - mT] which using eq.24 is w \^-mT . This is negative 
also. Thus both kinetic and potential energy fall giving rise to negative C p . 
If on the other hand we consider points a and b, point a has both a bigger 
volume and a bigger temperature thus 8m is positive. This would make both 
the kinetic and potential energy rise when one moves from b to a. This is 
illustrated in Table 1. The caloric curve of fig. 13 shows regions of negative 
C p . 
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Table 1 

Variation of energies per particle (MeV) with temperature (MeV) in the negative 
and positive compressibility zones, for p = 0.017 MeV fm~ s 
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Fig. 13. Caloric curve at constant pressure (p=0.017 MeV/m -3 ) in the canoni- 
cal model witfb4=200 and iV=200. The solid and dashed portions of the curve 
give -ve and +ve C p respectively. 



Let us consider the thermodynamic limit. This will be reached when the num- 
ber of composites near the boundaries of the freeze-out volume is negligible to 
the number of composites well inside. In this limit, intensive variables remain 
unchanged when extensive variables are changed by a constant factor. Thus if 
A is the total number of nucleons in the system and we change, at constant 
temperature, A — > A + aA, V — > V + aV the pressure p = yT must remain 
constant. This means, for constant T, when A — > A + aA, V — > V + aV, m 
must change to m — > m + am. Now for compressibility, A stays at A, but V 
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to V + aV thus m must change to less than m + am. Then the pressure will 
fall when V is increased, i.e., regions of negative compressibility disappear. 

It would be nice to demonstrate this feature directly by doing canonical cal- 
culations for larger and larger systems. The area over which negative com- 
pressibility appears does drop as larger and larger systems are used but the 
convergence is slow. Instead we will use the grand canonical ensemble to get 
to the A = oo limit. For a given density we solve eq.18, setting once k m = 200 
and k m = 2000, the other time. This means, in the first case, the largest 
composite has 200 nucleons and in the second case, the largest composite 
has 2000 nucleons. The temperature is chosen to be 6 MeV. Eq.18 has no 
reference to either A or V (only their ratio), the implication being for the 
grand canonical model to be good each factor is oo or very large. Pressure in 
the grand canonical ensemble is given by p = (T/V)lnQ gran d which leads to 
P = Efc=iexp(/c/i/3)^. 

Fig. 14 compares canonical calculation with A=200 and k m = N = 200 with 
A = oo and iV=200. We see in the low density (the gas phase) the two diagrams 
coincide. The rise of pressure with density is quite rapid and linear. After the 
two diagrams separate, the rise of pressure with density in the grand canonical 
ensemble slows down considerably but there is no region of mechanical insta- 
bility although the canonical calculation with 200 particles has a region of 
instability. In the grand canonical result which represents the thermodynamic 
extrapolation, we have not reached the classic liquid-gas coexistence limit 
where there would not be any rise of pressure at all (such as in Maxwell's con- 
struction). We think the reason is this. The largest cluster is 200 which is not 
a big enough number. We now increase the largest cluster size to 2000. Now 
the coexistence region is very clear and there is unmistakable signature of first 
order phase transition. In the same figure we also show results of canonical 
calculation with A=2000 and iV=2000. The region of mechanical instability 
has gone down considerably but it has not disappeared showing that we have 
not reached the thermodynamic limit yet. 

Thermodynamics allows C p to become negative. The following well-known 
relation exists [2]: 



C p - C v = VT— 

K 

where a is the volume coefficient expansion and k is the isothermal compress- 
ibility given by 



V [ dT )p 
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For negative k, C p is less than Cy and can become negative. 
Using the equality (|^) p = — (fr) T (irV we can a ^ so wr ^ e 



C n C\ 



dp dV 



(25) 



This shows that C p can drop below Cy if the isobaric volume coefficient of 
expansion becomes negative which is the case in some regions of fig. 11. 
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Fig. 14. EOS at T=6 MeV in the two models. For the left panel, the largest 
cluster has iV=200 and for the right panel iV=2000. For the canonical calcula- 
tion, the left panel has ^4=200 and the right panel has v4=2000. For the grand 
canonical calculation A = oo. 



We leave now general considerations of phase transitions, specific heat, caloric 
curves etc. and explore the predictive powers of the canonical thermodynamic 
model in producing detailed data in heavy ion reactions. Specifically we will 
investigate how effective the canonical thermodynamic model is in predicting 
isotopic yields in some specific reactions. For this we need to go beyond the 
production of hot fragments that the canonical thermodynamics will give. To 
obtain yields of specific final products, we need to investigate how fragments 
at non-zero temperatures will decay. The next sections address this issue. 
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7 Corrections for secondary decay 

The statistical multifragmentation model described above calculates the prop- 
erties of the collision averaged system that can be approximated by an equi- 
librium ensemble. Ideally, one would like to measure the properties of excited 
primary fragments after emission in order to extract information about the 
collisions and compared directly with the equilibrium predictions of the model 
described in this report. However, the time scale of a nuclear reaction (1CT 20 
s) is much shorter than the time scale for particle detection (1CT 9 s) . Before 
reaching the detectors, most particles decay to stable isotopes in their ground 
states. Thus before any model simulations can be compared to experimental 
data, it is indispensable to have a model that simulates sequential decays. This 
turns out to be not a simple task. 

In this section, we follow the techniques of refs. [57,58] to calculate the sec- 
ondary decay. We identify some issues that can be accurately addressed and 
others that are less controlled and may contribute uncertainties that influence 
the final results. Later, we calculate the secondary decay of excited nuclei 
predicted by the statistical multi-fragmentation model and compare the final 
ground state yields to recent measurements. 



7.1 Levels and level densities 

To calculate the secondary decay corrections, one must specify both the high 
lying states that are mainly populated at freeze-out and the lower lying states 
whose populations increase as these excited nuclei decay towards the ground 
state nuclei that are experimentally measured. In the previous sections, the 
discussion was centered about the states that are populated at freeze-out. 
While in principle all nuclear states may be involved at freeze-out, the vast 
majority of fragments are excited to the particle unbound continuum. 

The level densities in the unbound continuum influence the overall yield of un- 
bound nuclei at freeze-out as well as the sequence and the number of particle 
decays. In principle, interactions between fragments and their surroundings 
modify the states and their excitation energies. The vanishing of the surface 
tension cr(T) in the free energy expression at the critical temperature T c = 18 
MeV reflects such considerations. Few experimental constraints on continuum 
level densities exist, however, even when the nuclei are isolated. Thus, un- 
certainties in the continuum level densities introduce uncertainties into the 
calculated results. 

Following ref. [57], we represent the continuum level densities corresponding 
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to the internal free energies in Eqs. 7 and 13 by the expression: 



P smm{E\ J) = p S MM(E*)f(J, a). (26) 



4o()0 , 



where Psmm(E*) = p FG (E*) e -b SMM ^ SMM E^/^ bsMM = o.07A" 1 - 82 ( 1 " 
asMM = ^ + §co4j^-, J is the spin, E* is the excitation energy and A is the 
mass of the fragment. For light and medium mass nuclei, clsmm ~ A/%. Here, 



rf) ° v^f/' eIP ( 2 ^»' P )- (27) 

(27 + l )M p[-(J + 1/2) Wl 
2o" z 



a 2 ~ 0.0888y A • i?*/8)A 2//3 , (29) 

and -E* and Z are the excitation energy and charge of the fragment. For further 
details, we refer the reader to ref. [57]. 

In contrast to the continuum level densities, the discrete level densities need 
no corrections for the influence of interactions because these levels become 
important only much later in the decay after the fragments have decoupled 
from their surroundings. For this purpose, we use the spectroscopic informa- 
tion of isolated nuclei with Z < 12 where the information is available. For 
12 < Z < 15, low-lying states are not well identified experimentally and a 
continuum approximation to the discrete level density [59] was used. For all 
fragments with Z < 15 and excitation energies between the domains of discrete 
and continuum level densities, the level densities were smoothly interpolated 
[57]. 

Where the experimental information for nuclei with Z < 15 is incomplete, 
values for the spin, isospin, and parity were chosen randomly in the decay 
calculations as follows: spins of 0-4 (1/2-9/2) were assumed with equal proba- 
bility for even-A (odd- A) nuclei, parities were assumed to be odd or even with 
equal probability, and isospins were assumed to be the same as the isospin of 
the ground state. This simple assumption turns out to be sufficient since most 
of spectroscopic information is known for these low-lying states. 

For excitation energies where little or no structure information exists, levels 
were assumed to be specified by the relevant level density expression. Groups 
of levels were binned together in discrete excitation energy intervals of 1 MeV 
for E* < 15 MeV, 2 MeV for 15 < E* < 30 MeV, and 3 MeV for E* > 
30 MeV to reduce the computer memory requirements. The results of the 
calculations do not appear to be sensitive to these binning widths. A cutoff 
energy of E* uto j f /A = 5 MeV was introduced corresponding to a mean lifetime 
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of the continuum states at the cutoff energy about 125 fm/c. Where unknown, 
parities of these states were chosen to be positive and negative with equal 
probability and isospins were taken to be equal to the isospin of the ground 
state of the same nucleus. In this fashion, a table of states for nuclei with 
Z < 15 was constructed. 



7.2 Sequential decay algorithm 



Before sequential decay starts, hot fragments with Z < 15 were populated 
over the sampled levels in the prepared table according to the temperature. 
For the ith level of a given nucleus (A,Z) with its energy E* and spin Jj, the 
initial population is, 

Y, = Y (A,Z) P/.+ l)e*P( 7 ^M^J.) 

W l Y.,(2J i + l)ex P (-Et/T)p(Et,J,) 1 ' 



where Yq is the primary yield summed over all states of nucleus (A,Z) and T 
is the temperature associated with the intrinsic excitation of the fragmenting 
system at breakup. 

Finally all the fragments will decay sequentially through various excited states 
of lighter nuclei down to the ground states of the daughter decay products. The 
decay of fragments with Z > 15 was calculated according to the fission model 
of ref. [60]. The subsequent decay of excited fission fragments with Z < 15 
was calculated according to the Hauser-Feshbach algorithm described here. In 
this algorithm, eight decay branches of n, 2n, p, 2p, d, t, 3 He and alpha were 
considered for the particle unstable decays of nuclei with Z< 15. The decays 
of particle stable excited states via gamma rays were also taken into account 
for the sequential decay process and for the calculation of the final ground 
state yields. If known, tabulated branching ratios were used to describe the 
decay of particle unstable states. Where such information was not available, 
the branching ratios were calculated from the Hauser-Feshbach formula [61], 

^ = (31) 



where 



Gd — (Idleldzle^lplpz) 2 

\J d +Je\ \J P +J\ 1+ /_ 1 y 
x £ £ ! + We( l) Ti{E) (32) 

J=\J d -Je\ l=\J p ~J\ 
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for a given decay channel d (or a given state of the daughter fragment). J p , 
Jd, and J e are the spins of the parent, daughter and emitted nuclei; J and / 
are the spin and orbital angular momentum of the decay channel; Ti(E) is the 
transmission coefficient for the lih partial wave. The factor [l + 7r p 7r d 7r e (-l)']/2 
enforces parity conservation and depends on the parities it = ±1 of the parent, 
daughter and emitted nuclei. The Clebsch- Gordon coefficient involving I p , Id, 
and I e , the isospins of the parent, daughter and emitted nuclei, likewise allows 
one to take isospin conservation into account. 

For decays from empirical discrete states and I < 20, the transmission coef- 
ficients were interpolated from a set of calculated optical model transmission 
coefficients; otherwise a parameterization described in Ref. [59] was applied. 



8 Comparisons to data 

Even though the structure of the low-lying states of the fragments plays little 
role in properties of the hot system, these structure effects become critical 
when the fragments cool later by secondary decay. In the sequential decay 
algorithm described in the last section, in addition to more sophisticated level 
densities, empirical binding energies of the known nuclei are incorporated. 
Where the empirical masses are lacking, an improved mass formula [16,57] 
is employed. To be self-consistent, the same masses and level densities are 
used both in the thermodynamic model which produces the excited primary 
fragments and in the subsequent sequential decay. This self-consistency re- 
quirement appears to be necessary [16] for some observables. The resulting 
code which combines the thermodynamic model with the sequential decay al- 
gorithm is referred to as ISMM for improved Statistical Multifragmentation 
Model in the following sections of the report. 

To illustrate the capabilities of the thermodynamic model, we calculate final 
ground state elemental and isotopic yields for systems with A = 168 and 
Z = 75 and A =186 and Z =75 at T = A.7MeV, corresponding to E*/A « 5 
MeV. In all the following calculations, the freeze-out density is taken to be 
1/6 of the saturation density. These two systems were chosen because they 
have the same proton fractions as the combined systems formed in central 
112 Sn+ 112 Sn and 124 Sn+ 124 Sn collisions, respectively. However, the overall size 
and excitation energy of these systems have been reduced below that of the 
corresponding compound nuclei to reflect the loss of particles and excitation 
energy to pre-equilibrium emission prior to the multi-fragment breakup. These 
parameters have not been adjusted to obtain a best fit of the data. 

In the following, we illustrate the capability of this thermodynamic model 
to describe experimental charge, mass and isotopic yield distributions. We 
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also compare experimental and calculated observables, such as the isotopic 
temperature and the isoscaling parameters, which are constructed from these 
yields. 



8.1 Charge and Mass Distributions 



Calculations of the mass distribution for excited primary fragments are shown 
in Fig. 15 for a system with Aq = 168 and Zq = 75 at T = 4.7 MeV. The 
distributions of the primary fragments directly obtained from the thermody- 
namic model are shown as dashed lines with open points while the solid line 
with solid points represent the distributions of the final fragments after se- 
quential decays. Certain differences between primary and final spectra can 
be expected. Heavier fragments formed in the multifragment stages decay to 
smaller fragments, shifting the distribution to lower masses. In addition, the 
decay produces a large increase in the hydrogen and helium particles, because 
these are the main products of the decay of the heavy fragments. 
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Fig. 15. Predicted mass distributions from the multifragment at ion of a source 
nucleus with the mass number 168 and charge number 75. The open circles 
are primary yields and the closed circles are yields after secondary decay. 

The differential multiplicities dM/dVt for various masses with A < 20 are 
plotted in an expanded scale in Fig. 16 for both the A = 168 and A = 186 
systems. For comparisons, experimental data obtained by averaging over 70° < 
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Ocm < 110° for central 112 Sn+ 112 Sn and 124 Sn+ 124 Sn collisions at E/A=50 MeV 
[62] are plotted as open and solid points in the left and right panels, respec- 
tively. The calculations reproduce many features of the mass distribution. 
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Fig. 16. Predicted mass distributions (A < 20) from the multifragmentation 
of a source nucleus with A =168 and Z =75 (left panel) and A =186 and 
Zq=75. The dashed lines are the predicted primary yields and the solid lines 
are predicted yields after secondary decay. For comparison, data from multi- 
fragmentation of central collisions of 112 Sn+ 112 Sn are shown as open symbols 
(left panel) and closed circles for 124 Sn+ 124 Sn reaction (right panel) [62]. 

The relative normalization of the calculation can be increased by increasing 
the size of the source or by making its angular distribution sideways peaked. 
The slope of the mass distribution can be made more steep by increasing 
the source temperature. There are indications that the experimental angular 
distributions are not isotropic and that pre-equilibrium emission mechanisms 
may contribute to the yields of the lighter fragments. Accordingly, we do not 
fit the calculations to the experimental data in this article, but defer such 
detailed analyses until more experimental data that can constrain such effects 
become available. 

The charge distributions exhibit similar behavior as the mass distributions. 
For completeness, we include the charge distributions for the Aq = 168 and 
Zq = 75 and Aq = 186 and Zq = 75 in Figs. 17 and 18. The same conventions 
for the mass distribution figures (Fig 15 and 16) are used. 
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Fig. 17. Predicted charge distributions from the multifragmentation of a source 
nucleus with v4o=168 and Zq=75. The open circles are primary fragment yields 
and the closed circles are after secondary decay. 
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Fig. 18. Predicted charge distributions [Z < 8) from the multifragmentation 
of source nuclei with A =168 and Z Q =75 (left panel) and A =186 and Z =75. 
The open and solid points are data from Ref.[62]. See Fig. 16 for explanations 
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of symbols used. 



In the break up calculations, the odd-even effects are evident. These occur 
because pairing and shell effects are not completely washed out in our level 
density expressions at a temperature of T = 4.7 MeV. As the secondary decay 
washes out such structures, these odd-even effects in the primary distribution 
have little or no effect on the final fragment distribution. 



8.2 Isotopic distributions 
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Fig. 19. Isotope distributions for Carbon and Oxygen fragments. The dashed 
and solid lines correspond to the predicted primary and final yields respec- 
tively. The open and solid points are from Ref. [62]. 

In Fig. 19, the isotopic distributions for carbon and oxygen isotopes are plotted 
for the two sources. Using the same convention as before, the dashed lines 
correspond to the distributions of the primary fragments while the solid lines 
represent the final distributions after sequential decay. As expected, the more 
neutron-rich system with N /Z Q = 1.48 preferentially produces more neutron- 
rich isotopes than the neutron deficient system with N /Z Q = 1.24. In all cases, 
the primary distributions are much wider and more neutron-rich than the final 
distributions. The experimental isotope distributions (data points) agree more 
with the final results obtained after secondary decay than with the primary 
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distributions. Nonetheless, the widths of the experimental distributions exceed 
those of the final distributions and are more neutron-rich. This suggests that 
the predicted corrections for secondary decay may be somewhat too large. 
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Fig. 20. The mean neutron to proton ratios as a function of the charge of the 
emitted fragment Z for the two systems. The left and right panels correspond 
to the calculated results from the primary and final fragments. 

The mean neutron to proton ratios < N/Z > for each element provides another 
observable with sensitivity to the isospin asymmetry dynamics of the reaction. 
The dependence of the calculated primary values on the < N/Z > of the total 
system is much stronger than that of the final values. This can be seen in 
Fig. 20 where the primary (left panel) and final (right panel) < N/Z > values 
are compared for the two systems. The differences of the primary values for 
< N/Z > of the two systems are large, reflecting the large difference in the 
initial isospin asymmetry of the two systems. The largest values for < N/Z > 
occur for Z ~ 8,20, etc., values corresponding to nuclei where one can have 
either closed proton or neutron shells. Such nuclei can remain comparatively 
well bound even for large value of N/Z. Both of these enhancement and the 
difference between the < N/Z > values for the two systems are diminished 
in the final distributions, which are both narrower and located closer to the 
valley of beta stability. 

Fig. 21 shows measured and calculated primary and final values for < N/Z > 
as functions of the element number Z. The left and right hand panels provide 
the < N/Z > values for the neutron-deficient and neutron- rich systems, re- 
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spectively. The calculated final distributions reproduce the measured values 
well. It is rather curious that the experimental < N/Z > values exhibit the 
odd and even effects as a function Z. Such staggering is much less obvious in 
the neutron rich system. For reference, the < N/Z > for the abundances of 
naturally occurred isotopes are plotted as stars in both panels of the figure. 
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Fig. 21. The mean neutron to proton ratios as a function of the charge of 
the emitted fragment Z for the neutron deficient (left panel) and neutron rich 
(right panel) systems. For comparison, data from the multifragmentation of 
central collisions of 112 Sn+ 112 Sn are shown as open symbols (left panel) and 
as closed circles for 124 Sn+ 124 Sn reaction (right panel) [62]. For reference, the 
mean N/Z ratios from naturally occurring isotopes are shown as stars. 



8. 3 Isoscaling 

The dependence of the isotopic distributions on the N /Z of the colliding sys- 
tem can be more sensitively explored by the use of isotopic ratios [62,63,64,65]. 
In particular, the ratio, R,2i(N,Z) = Y 2 (N, Z)/Y 1 (N, Z), of yields from two 
different reactions, labelled here as 1 and 2, has been shown to exhibit an 
exponential relationship as a function of the isotope neutron number N, and 
proton number, Z [62,63,64,65,66,67,68,69,70,71,72,73,74]. 

R 2 i(N,Z) = C ■ exp(aN + (3Z) (33) 



N /Z =1.24 


i i i i 
/- X N /Z = 1.48 


★ 






/ *\ 




/ Jyju; 






- // ★ ★ 


/ ★ ★ 


\ primary 




final 


- * valley of stability - 


- / o 112 Sn+ 112 Sn 

, , , i , , , , i , , , , i , , , , i 


- 124 c . 124 c 

• Sn+ Sn 

, , , i , , , , i , , , , i , , , , i 



41 



where C is a normalization factor and a and (3 are the isoscaling parameters. 

Calculations with a variety of different statistical models show that the isoscal- 
ing relationship is strictly obeyed by the primary fragments in these models 
[64,66,71]. Surprisingly the isoscaling relationship is also obeyed by fragments 
produced in dynamical models such as the asymmetrized molecular dynamical 
model [70]. In all cases, the isoscaling parameters are related to the isospin 
asymmetry of the collisions and to the form of symmetry energy or, equiva- 
lently, asymmetry term of the EOS chosen in the model [64,66,70,71,75]. 

Neglecting for simplicity the Coulomb interactions between fragments and 
environment, the exponential dependence of the isoscaling relationship can 
be easily understood from the expression for the yields for a fragment with 
neutron and proton numbers N and Z within the grand canonical limit of the 
present equilibrium model [76]: 

Yi(N, Z) = y / 3/2 y (T) exp [(Z fi p , + N» nti + B N)Z )/T] . (34) 



Here, qN,z{Ti) represents the internal partition function of the fragment, Vi 
the free volume of the system, At = ^j2irh 2 /mTi, m the nucleon mass and ji pi 
(// n> j)the chemical potential associated with free protons (neutrons) for the 
i th reaction which produces a system at temperature Tj. If the temperature 
in the two reactions are expected to be the same (as in the Sn reactions 
described here), the chemical potentials ji pi and ii n>i contain the only reaction 
dependent factors in this exponential. In this limit, a = \\i n ,i ~ fi>n,i]/T and 

The symbols in Fig. 22 represent the isotopic ratios calculated by the canon- 
ical thermodynamic model described in this review. In Figs. 22 and 23, the 
following convention is adopted. We choose closed symbols and solid lines for 
even Z and open symbols and dashed lines for odd Z starting with Z=l for 
leftmost line. The lines are best fits of the calculated R21 ratios to Eq. (33); 
the lines are essentially linear and parallel on this semi-log plot consistent with 
a single constant isoscaling parameter a pr i mary = 0.50. The spacing between 
these lines corresponds to the increase in i? 2 i for unit increases in Z; the ob- 
served equal spacing is consistent with a single constant isoscaling parameter 

ft primary 0.64. 

For comparison to the data, we only examine the isotope ratios where there 
are data with sufficient statistics. The symbols in the bottom panel of Fig. 23 
represent the predicted isotopic ratios after sequential decays. The lines are 
nearly parallel to the lines in Fig. 22 on average and the isoscaling parameters 
Oifinai = 0.46 and (3fi na i = —0.52 are comparable to the primary values. In 
detail especially when the isotopes away from the valley of stable nuclei are 
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considered, the trends are not as clearly consistent with the isoscaling law as 
are the trends of the primary distribution. The larger change in the (3 values 
may arise from the approximation of the Coulomb interaction used in the 
model. In the top panel, the data are shown as symbols. The experimental 
isoscaling parameters are a data = 0.36 and fldata = —0.42. The slopes from the 
calculations are flatter suggesting that the temperature of 4.7 MeV used as 
the input parameter in the model may be too low. However, if the temperature 
is increased so that the isoscaling predictions agree with the data, the other 
observable such as the mass and charge distributions as well as the isotope 
distributions may no longer agree. As stressed earlier, the current work is 
not to use the optimized set of model parameters but rather to compare the 
trends of data with the the model calculations. More constraints and study 
are needed to optimize the agreement with data. 




Fig. 22. Predicted yield ratios, R 2 i(N, Z) = Y 2 (N, Z)/Y 1 (N, Z) from primary 
fragments for the two systems studied in this work. The lines are best fit to 
the symbols according to Eq. (33). Different lines correspond to Z=l to 8 
starting with the leftmost line with three points being Z=l. 
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Fig. 23. Top Panel: Experimental isoscaling behavior exhibited by the central 
112 Sn+ 112 Sn and 124 Sn+ 124 Sn collisions. The data are the nuclide yield ratios, 
B,2i(N, Z) from the two reactions plotted as a function of N. The isotopes of 
different elements lie along different lines. The solid and dashed lines repre- 
sent the best fit to Eq. (33). Bottom Panel: Predicted yield ratios, R2i(N, Z) 
obtained from the final yields for the two systems studied in this work. The 
symbols and lines have the same convention as the data used in the top panel 
and Fig. 22. 



8.4 Isotopic temperatures 



Starting from the grand canonical expression for the yields (Eq. 34), it is 
also possible to construct a double ratio that minimizes the sensitivity to the 
isospin asymmetry while maximizing the sensitivity to the temperature. By 
doing so, one can construct an isotopic thermometer, whereby the temperature 
is extracted from a set of four isotopes produced in multifragment breakups 
as follows [77] , 

rp Aff 

Ilso ~HaR) [6b) 



where 
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AB 



R 



Y(A 2 ,Z 2 )/Y(A 2 + 1.Z 2 ) : 
B(A 1 ,Z 1 ) - B(A 1 + l,Zi) 



(36) 



-B(A 2 ,Z 2 ) + B(A 2 + 1,Z 2 ), 



(37) 



and 



(2J Z2 , A2 + 1) (2J ZuAl+1 + 1) f A 2 (A, + 1) 
(2J ZliAl + 1) (2J Z2iA2+1 + 1) [A, (A 2 + 1) 



n3/2 



a 



(38) 



In this ratio derived from Eq. (36) for the ground state yields, Y(A, Z) is the 
yield of a given fragment with mass A and charge Z; B{A, Z) is the binding 
energy of this fragment; and Jz,a is the ground state spin of the nucleus. In 
the context of the grand canonical ensemble, Eq. (35) has been regarded as an 
effective or "apparent" temperature that may differ somewhat from the true 
freezeout temperature T due to the influence of secondary decay and other 
cooling mechanisms. 

The influence of secondary decay on the isotopic temperatures can be clearly 
observed because it leads to variations in the values for the temperature that 
depend on the isotopes used to construct the ratio. The variations are uni- 
versal, observed in many different reaction systems and thus can be used to 
assess the effectiveness of sequential decay models. One origin of these varia- 
tions is the feeding from higher lying particle bound states. Such effects can 
be modeled by changing the value for the statistical factor "a" and making it 
temperature dependent. This and additional feeding from the decay of heav- 
ier particle unbound nuclei can be modeled by the secondary decay formalism 
described in the previous section. 

To illustrate the influence of secondary decay on isotope temperature mea- 
surements, measured and calculated final temperatures have been extracted 
from double ratios of Z=2-8 fragments and plotted in Fig. 24. To reduce the 
influence of secondary decay, we include only isotope thermometers with large 
values for AB in this figure. This requirement restricts comparisons to three 
types of thermometers: a.) T iso ( 3 > 4 He) with Z 2 =2, A 2 =3, b.) T iso ( u ' 12 C) with 
Z 2 =6, A 2 =U, and c.) T iso ( 15 - 16 0) with Z 2 =S, A 2 =15. We note that the ther- 
mometer (a) involves the light particle pair 3 ' 4 He while thermometers (b) and 
(c) concern only intermediate mass fragments with Z = 3 — 8. The solid lines 
show corresponding ISMM predictions for these three types of thermometers 
as a function of A\. 

Similarities in the variations of the calculated and measured temperatures 
allow insight into their origin. Each panel of Fig. 24 corresponds to fixed values 
of Z 2 and A 2 ; the observed variations in T iso are therefore correlated with Z\ 
and A x . The highest values for T iso involve 10 Be (Z 1= 4, A 1 = 10) and 18 (Z 1= 8, 
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A\ = 18). The calculations attribute this increase to enhancements in the yields 
of these nuclei due to 7-ray feeding from their many low lying particle bound 
states [80]. Other thermometers in Fig. 24 provide temperature values that 
are significantly lower than those involving 10 Be and 18 O. Most thermometers 
are significantly lower than the primary temperature of 4.7 MeV, depicted by 
the horizontal dashed line in the three panels. 
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Fig. 24. Isotopic temperatures extracted from three types of thermometers. 
Experimental data are shown as symbols. The lines are calculations. For ref- 
erence, the input primary temperature of 4.7 MeV is shown as the horizontal 
dashed lines. (See text for details on the dot-dashed line in the left panel.) 

Both calculated and measured values display a Z or A dependence in T iso . 
Calculated values for Tj SO ( 15,16 0) are about 0.5 MeV lower than those for 
7iso( n ' 12 C), which are about 0.2 MeV lower than T iso ( 3,4 ife). There is also 
a trend for isotopic temperature values to decrease as a function of A\. The 
calculated decrease of T iso with A\ and A 2 reflect the increasing importance of 
multi-step secondary decay contributions to the yields of these heavier nuclei. 
Such multi-step decays make the system appear cooler because the final ground 
state nuclei originate from the decay of an ensemble of unstable nuclei that 
are less excited than the original ensemble. 

We note that the experimental Ti SO ( 3,4: He) temperatures (solid symbols in 
the left panel) are systematically higher than the corresponding ISMM values 
(solid line). As these thermometers derive their sensitivity to the temperature 
from the large binding energy difference between 3 He and 4 He, the difficulty in 




46 



reproducing these quantities may arise if there are significant pre-equilibrium 
production mechanisms for light particles such as 3 He [80]. To illustrate this 
effect, we assumed that 2/3 of the measured 3 He yield is of a non-thermal 
origin. This increases the 3 He yield by a factor of three; calculations including 
this pre-equilibrium enhancement are shown as the dot dashed line in the left 
panel. The success of this resolution of the discrepancies between Ti SO ( 3,4 He) 
and Tj SO ( 11,12 C) suggests that it may be necessary to make careful estimations 
of the contributions from pre-equilibrium emission before isotope temperature 
measurements involving T iso ( 3 ' 4 He) will be fully accurate. 



9 Summary 

The canonical version of the thermodynamic model has helped clarify many 
aspects of intermediate energy heavy ion collisions. The obvious advantage 
is that, as opposed to the grand canonical model, it has an exact number of 
particles. The predictions of the grand canonical model (which really applies 
to very large systems) can differ very significantly from those of the canoni- 
cal model specially in the intermediate energy regime. The canonical model 
helps us to understand the order of phase transition, the caloric curve and 
the possibility of negative specific heat. The model gives quantitative fits to 
experimental data on isotopic yields and the phenomenon of isoscaling, now 
so well established in intermediate energy heavy ion collisions. The virtue of 
the model is also its simplicity. Most of the calculations reported in this work 
can be carried out quite easily. 
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A Equilibrium, reactions and reaction rate time scales 

A basic assumption of statistical models is that equilibrium is reached in the 
time scale of the reaction. For fragment or composite particle distributions a 
complex set of reactions takes places [81,82]. The processes involved in the 
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collision of heavy ions can be modeled in a manner that is similar to nucle- 
osynthesis in a dense, heated and evolving system such as in the expansion 
of the early universe and in supernovae explosions. The starting point of such 
a description is then a dense and heated system of neutrons and protons 
which combine through a set of reactions to make the composite nuclei from 
the lightest nuclei such as deuterons, alpha particles,.. , all the way up to 
much heavier and complex nuclei. By way of illustration and also for con- 
trast, the nucleosynthesis in the early universe occurs through a set of two 
body reactions with the first element of the chain being an electromagnetic 
radiative capture of a neutron plus proton to a deuteron with an emitted pho- 
ton carrying away the excess energy. After this first electromagnetic process, 
light elements are produced by a sequential set of two body reactions such as 
d + d — > He 3 + n, d + d — > t+p, t + d — > He 4 + n,... Nuclei up to Li are believed 
to be produced at their equilibrium concentration in big bang nucleosynthesis 
models. The abundance of heavy elements comes from processes involved in 
supernovae. The study of these processes is the area of nuclear astrophysics 
and heavy ion collisions offer the opportunity to study similar processes and 
phenomena in the laboratory. 

In heavy ion collisions, electromagnetic processes are too slow over the time 
scale of the collision to produce the observed distribution of composites or 
produced particles. A typical time scale of the collision is 10~ 22 sec or 30 fm/c 
which is much shorter than any electromagnetic process time scale. Densities 
in heavy ion collisions can be high enough for a three body process to occur 
such asn + p + N^-d + N, where the nucleon N can carry away the excess 
energy. At very high energies, meson production processes occur, so that a d is 
formed in radiative pion emission of n + p. Heavier composite particles evolve 
through reactions such as those listed above. However, it should be noted that 
because of possible very high initial densities, multi body processes can occur 
besides two body processes even for composites heavier than the deuteron. 
These only enhance the approach to equilibrium. At RHIC energies, particle 
production becomes very important, and reactions leading to new particles 
have been studied [82,83]. 

As an example of a reaction rate approach consider the formation of a deuteron 
through the process p + n + N^d + N. The time evolution of the deuteron 
density pd can be obtained from an equation involving the proton density p p , 
neutron density p n , and nucleon density p^: 

^ L = [p P Pn[(—)e q -Pd}pN(a[N + d^n + P + N}xv) (A.l) 

at PnPp 

Here (-fy)eq is the equilibrium ratio of the densities of gPs to ra's and p's 
and is a function of temperature. The () term in A.l involves the product 
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of the breakup cross section of deuterons induced by nucleons and v, which 
is the relative velocity of the N and d pair. This product is averaged over 
the velocity distribution of the pair. In obtaining the expression in A.l we 
used detailed balance which relates the forward rate for the formation process 
p + n + N — > d + N to the backward rate of the break up or absorption 
process d + N — > p + n + N. Equilibrium is reached when the forward rate 
is equal to the absorption rate. Initially, the deuteron density is being built 
up by forward processes which involves the product of proton, neutron and 
nucleon densities, but later in this time evolution deuterons will start to be 
absorbed by backward processes which involve the newly formed deuterons and 
the existing nucleons. Once equilibrium is reached these underlying processes 
vanish in the description of the deuteron density, which is now described by 
phase space factors with temperature and volume playing a dominant role. 
Large volumes reduce composites since nucleons are less likely to be near 
each other to combine and high temperatures increase break up probabilities. 
Binding energy terms appear as Boltzmann factors and enhance composite 
densities. 

We can question whether rates are fast enough to produce equilibrium dis- 
tributions. To answer this question we consider the following simplified ex- 
pression for a reaction rate: pn x cr x v . For p^ we take nuclear matter 
density or 0.15 nucleons/fm 3 . Typically temperatures are 10's of MeV for 
medium energy collisions and a temperature of lOMeV has a kinetic energy 
of 15 MeV = (l/2)m(t>/c) 2 . For v/c = 1/5, a cross section lfm 2 will 
have a rate 10 22 /s. The reciprocal of this rate is the reaction rate time scale 
which is 10 _22 s. Thus, a cross section of lfm 2 will have a reaction rate time 
scale that is equal to the characteristic time scale of the collision. Under these 
circumstances equilibrium will be reached. 

Next, consider the prototype two body reaction A + B — > C + D. The rate 
of growth of the density of C can be related to the chemical activity A = 
Pa + Pb — Pc — (J>d, where pa is the chemical potential of A, etc. Specifically, 
the time evolution of the density of C is 



= p APB (a[A + B ^C + D]xv}(l- expl-A/T]) (A.2) 

At equilibrium pa+Pb = Pc+^d- Thus, the factor (1— exp[— A/T]) — > 0. Near 
equilibrium A «T and (1 — exp[— A/T]) — > A/T. In this limit the reaction 
rate eq.(A.2) is linear in the chemical activity A. Such linear connections are 
known as Onsager relations where the chemical activity acts as a generalized 
force, X, and the left hand side of eq.(A.2) is interpreted as a generalized 
velocity J. Then J = LX, where L is the proportionality constant between J 
and X. Far from equilibrium, this linear relation is no longer valid since A is, 
in general, not small compared to T. 
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As a final consideration in discussing reaction rates we note that if the equi- 
librium concentration of the particle of interest is small, then the reaction rate 
constant is somewhat more complicated than the simplified expression used 
above. To illustrate this situation we mention the case of pion production. For 
example, for the reaction N + N^N + N + ir, the rate equation for the pion 
density is: 



Here,(p„-) e g is the equilibrium pion density which depends on temperature. 
This rate equation can be solved to give p n (t) = {p-w)eq * (1 — exp[— A x t\). 
The rate constant is 



The result of eq.(A.4) differs from the simplified reaction rate used above by 
an important factor Pn /{.P-^eq- This factor can be very large when the equilib- 
rium density of pions is small compared to the nucleon density. It was one of 
the reasons why the results of [82] led to the conclusions that pions would be 
in chemical equilibrium, a result which differed from a previous result in [84]. 
While low equilibrium concentration can enhance reaction rate constants and 
reduce equilibration time scales, some examples of other enhancement factors 
are the presence of two or more channels to the final state, the presence of 
secondary processes, high densities which allow multiparticle production pro- 
cesses above the two body type just considered. For example, the time scale 
for kaon production is considerably reduced through pion induced reactions, 
where the pions are copiously produced in the initial nucleon-nucleon collisions 
as first noted in [82]. 



B Antisymmetry and all that 

Our whole discussion started from eq.(2) in section 2 which then led to eq.(5), 
the recursive formula. Eq.(2) is not quantum mechanical. The partition func- 
tion of rii particles takes this simple form only under situation of low density 
and high temperature. We argue here that the approximation is quite good 
for intermediate energy heavy ion collisions. 

We start with qualitative arguments. The volumes used are about three times 
or more of the normal volume. At low temperature (~ 4 MeV) where one might 
imagine the approximation to fail, it survives because many composites appear 



d(p n )/dt = [p 2 N - (p 2 N *p 7r /(p 7r ) eg ] x (av). 



(A.3) 



A = (av) x p 2 N /{pir) 



eq 



(A.4) 
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thus there is not enough of any particular species to make (anti)symmetrisation 
an important issue. At much higher temperature the number of protons and 
neutrons increase but as is well-known the n\ correction takes the approximate 
partition function towards the proper one at high temperature. In a hypothet- 
ical world, the problem could get very difficult. Such a scenario would arise 
if the physics was such that at low temperature we only had neutrons and 
protons and no composites. An even worse situation would be if we had only 
neutrons (or protons). With these preliminaries let us proceed to estimate 
quantitatively the errors involved in actual cases that one might encounter in 
intermediate energy heavy ion collisions. 

The recursive relation eq. (5) is not limited to the approximation of eq.(2). 
It is shown in [18] that by regarding the grand partition function (in our 
case this grand partition function incorporates correct (anti) symmetry among 
particles) as the generating function of the canonical partition function one 
derives a relation like eq.(5) 



1 N 

QN((3) = -J2 kx kQN-kW) (B.i) 

fc=l 

where Xk is not a one-particle partition function but is to be obtained from 
an expansion of the grand partition function. We illustrate this with first the 
example of only protons filling up orbitals i, j, k, ... in a box. Now 



lnQ gr (P, /i) = ]T ln{\ + e^" fe ) (B.2) 

i 



EE 



-e 



3 



The coefficient of e 13 ^ is x^ which then gives x& = ^— ^ — J2i e~ kl3ei When this 
expression for x^ is used in eq. B.I it generates the correct partition func- 
tion. Orbitals are given occupancies greater than one and then eliminated by 
subtraction. This can lead to severe round-off errors when applied to degen- 
erate Fermi systems but will not affect the application we envisage here. The 
number of of protons is given by 



xxQz-i 2x 2 Qz^ 2 Zx z Qq 

z = [ ^r + ^^ + -^r )lz (B - 3) 

The value of Q is 1. 

Anticipating generalisation we will call Xk in the above case yf^. The subscript 
1, means it is a "composite" with one proton and no neutron. The superscript 
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k means it is obtained from the k—th term in the expansion; yfj^ will contribute 
to x kfi . 

If instead we had a boson, deuterons for example, we would have 

ln[Q gr . can ((3, fip, fi n )] = ]T -ln(l - e^+^e-te) (B.4) 

i 

= ^ - e j ^ p+l3lXn ~^ (B.5) 

i j 3 

Thus in the case of deuterons y[ fc | (which would contribute to x^k) is given by 

We can treat an assembly of protons, neutrons, deuterons, tritons...etc. The 
recursive relation if the dissociating system has Z protons and N neutrons is 



Qz,n = y ix itj Qz-i,N~j (B.6) 

Z i=l,Z,j=0,N 

The average number of a composite with i\ protons and i 2 neutrons is given 
by 



^ 1^11,12 ^ Vii,i 2 Q Z— ii,N— ii /Qz,N + 2y\;\ i2 Qz -2ii,N~2i 2 /Qz,N + ••• (B.7) 

Unless one is in an extreme degenerate fermi system, one can evaluate the y 
factors by replacing sums with integration. For example, yf^ = — J2i e~ nf3ei 
where the sum is replaced by / e~ nf3t g(e)de = 2^(^p) 3 / 2 . Here V is the avail- 
able volume. We have included the proton spin degeneracy; m is the proton 
mass. For the deuteron, yf\ = \ \e~ k ^g{e)de. This is 3 x 2 3 / 2 ^(^f ) 3 / 2 ^ 
where Et, is the binding energy of the deuteron. It is clear how to compute 
contributions from other composites. 

We test the accuracy of the yields as calculated throughout the main text by 
comparing with a calculation where the complete theory of symmetrisation 
and antisymmetrisation is used. Subject only to the approximation that sum- 
mation over discrete states has been replaced by an integration over a density 
of states, the calculation is exact. The results are taken from [18]. We take 
the dissociating system to have Z=25 and N=25. The lowest temperature 
considered is 3 MeV (one might argue that at lower temperature a model of 
sequential decay is more appropriate). The highest temperature shown is 30 
MeV. We take a freeze-out volume in which the composites can move freely as 
three times the volume of a normal nucleus with 50 nucleons. Aside from neu- 
trons and protons we allow the possibility of composites. Excited states of the 
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Table B.l 



Comparision of claculations of average yields and E/A. By exact we mean a calcu- 
lation with proper symmetry. Sum over discrete orbitals in a box has been replaced 
by integration as is the usual practice. 



Calc 


P 


n 


d 


t 


3 He 


4 He 


Z > 12 


Temp. 

(MeV) 


E/A 
(MeV) 


approx 


0.307 


0.032 


0.050 


0.007 


0.054 


0.679 


0.945 


3 


-7.863 


exact 


0.306 


0.031 


0.051 


0.007 


0.053 


0.696 


0.945 


3 


-7.861 


approx 


1.174 


0.898 


1.177 


0.560 


0.641 


2.489 


0.051 


6 


-4.117 


exact 


1.117 


0.856 


1.195 


0.553 


0.638 


2.573 


0.050 


6 


-4.135 


approx 


4.127 


3.955 


4.812 


2.099 


2.052 


1.985 


0.000 


12 


4.401 


exact 


3.860 


3.696 


4.941 


2.090 


2.051 


2.021 


0.000 


12 


4.308 


approx 


10.937 


10.893 


7.664 


1.686 


1.650 


0.379 


0.000 


30 


28.914 


exact 


10.512 


10.468 


7.885 


1.732 


1.696 


0.395 


0.000 


30 


28.844 



composites were not allowed (they could have been included but the purpose 
of the exercise was to compare two models: calculations without the inclusion 
of excited states were sufficient to reach conclusions). Spins and binding ener- 
gies for deuteron, triton, 3 He and 4 He are taken from experiments. For higher 
mass composites the binding energy is taken from empirical mass formulas. 
For fermions, spin 1/2 was assumed and for bosons spin was assumed. For 
each Z we take N = Z — 1, Z, and Z+l. We present in the table average yields 
of protons, neutrons, tritons, 3 He, 4 He and the sum of yields of all nuclei with 
charges greater than 12. The temperature range of 3 to 6 MeV are of interest 
to many experiments. We also show results at 30 MeV. The approximation 
used in the main part of the text is seen to be quite good. 



C Applications to other areas 

While the main emphasis of this report is on the thermodynamic model for 
nuclear multifragmentation , the applications of the approach developed in 
section 2 to other areas will be mentioned in this appendix. In particular, 
many problems in statistical mechanics can be reformulated in terms of equa- 
tions 1-5 in that section. Each problem has a different choice for the factor 
uj that appears in these equations and a different interpretation of it within 
the general structure of those equations. We will now illustrate these remarks 
with some examples. 

Let us consider the following parallel between multifragmentation and per- 
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mutations, which appear when Fermi-Dirac and Bose-Einstein statistics are 
included into problems with identical particles. Any permutation can be bro- 
ken up into cycle classes and this cycle class decomposition is the basis for 
this parallel. A given permutation of A particles has a specific cycle class de- 
composition which specifies the number of cycles of length k. This number is 
similar to the number of clusters of size A; in a fragmentation. Moreover, the 
same type of sum rule holds as with clusters. That is, for any given permu- 
tation, the total A is equal to the sum of the cycle length times the number 
of cycles of that length in that specific permutation. The canonical partition 
function for non-interacting particles such as Fermi-Dirac or Bose-Einstein 
particles in a box or in a one body potential well such as a harmonic oscillator 
well has a form given by eq.2 in section 2 [85,86,87]. For identical particles 
in a box of volume V and a system at temperature T, the uj weight factor 
for a cycle of length k is that of eq.6 with the q k = 1 in that equation for 
Bose-Einstein particles and q k = (— for Fermi-Dirac particles. Once the 
canonical partition function is obtained from the recurrence relation of eq.5, 
the thermodynamic free energy F can be calculated and all other thermody- 
namic quantities also follow from F. For example, the pressure will have a 
form involving an expansion in density and quantum volume which gives the 
quantum corrections to the ideal gas law coming from the symmetrization or 
anti-symmetrization of the particles. Bose-Einstein particles in a laser trap 
which is taken as a harmonic oscillator well have also been studied using this 
approach [87]. Fermions in a well can also be studied as mentioned in [87] and 
an extended discussion can be found in [88]. Interactions can also be included 
along with quantum statistics as shown in [18]. Some further observations re- 
garding permutations are as follows. The result of eq.4 gives the mean number 
of cycles of length % in terms of the ratio of the two partition functions A — i 
and A, and the u factor for that length. Near the Bose-Einstein condensa- 
tion transition long cycle lengths start to appear and this manifestation of the 
transition is analogous to the appearance of large clusters around the liquid 
gas phase transition. The results of eq.4. give the probability of a particular 
permutation, specified by its n vector, being present. In RHIC collisions many 
pions are produced and the application of the methods in sec. 2 can also be 
given. For example Bose-Einstein effects associated with thermal pions have 
been studied in [89,90]. For thermal pions at temperature T in a volume V 
the cycle length u factor of eq.6 is given by, 




Here, m is the mass of the pion and K 2 is a MacDonald function. The u 
weight factor also appears in expressions concerning the mean number of pions, 
its fluctuations, and in higher moments of the pion probability distribution. 
Examples of these connections are: 
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(N 2 )-(N) 2 = J2k 2 iu k 
<(JV-<JV)) S )=E*V 



(C.l) 



The sums that appear in eq.C.l are over all fc's. Note that Poisson statistics has 
only unit cycles, or k — 1 only in the sums. Then (iV 2 ) — (iV) 2 = (N). The pres- 
ence of cycles of length 2 and higher cycles produces departures from Poisson 
statistics. An important observation related to Poisson statistics comes from 
the fact that coherent states have associated Poisson distributions. Moreover, 
departures from Poisson statistics are associated with chaotic emission pro- 
cesses. At high temperatures, Maxwell-Boltzmann statistics apply which leads 
to Poisson statistics in statistical models. The pion probability distribution for 
having N pions is the ratio of the canonical partition function for a system 
of size N divide by the grand canonical partition function. This probability 
was investigated in [89] for the case of 158 GeV Pb + Pb collisions where it 
is shown to have a Gaussian shape with a width that is about 10% larger 
than a Poisson distribution with the same mean number of pions. Many other 
models of pion and, in general, particle multiplicity distributions can be de- 
veloped in a similar manner by specifying another form for uj k . Once uj k is 
given, all quantities of interest follow. The importance of a phenomenological 
approach to multiparticle distributions, which is based on known distribu- 
tions from probability theory, is shown in [91,92,93,94,95]. Moreover, a wide 
range of physical processes can be accommodated using such an approach. A 
specific and frequently used distribution is the negative binomial distribution 
where uj k = x%-. The symbol x is the negative binomial parameter while t 
is another parameter that is important in fixing the mean number of pions 
and its variance: (N) = x-^ and (N 2 ) - (N) 2 = (N)(l + ((N)/x)). The 
generalized approach in [89,90] also includes several well known specific prob- 
ability distributions as special cases of a more general distribution. Here, we 
will just mention a few examples of various phenomena that can be found in 
[89,90] which are as follows: (1) Emission from systems with a variable signal 
to noise ratio, where the signal is related to a Poisson processes which may 
originate from a coherent state and a noise level given by a negative binomial 
distribution. (2) Field emission from Lorentzian line shapes and its connection 
to a Feynman- Wilson gas [96]. (3) Pion laser models [97,98] and the role of 
Bose-Einstein enhancement for a Poisson emitting source. (4) Multiparticle 
emission as a one dimensional random walk process along a jet axis. A reader 
interested in the application of the methods of sect. 2 to multiparticle mul- 
tiplicity distributions can find the details and several other individual cases 
in [89,90]. In a series of papers [85,86], Hegyi has considered many interest- 
ing aspects of multiparticle production and has also introduced a generalized 
distribution for its description. 
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Photon count distributions can also be developed using the approach of sec- 
tion 2. In fact, early models of pionic distributions [91] coming from nucleon- 
nucleon and nucleus-nucleus collisions were based on photon count distribu- 
tions [91]. The laser distribution of [91] is an example of a distribution which 
first appeared in quantum optics and was then subsequently taken over into 
the area of particle production. Thermal emission of photons have an ujk factor 
that can be obtained as the zero mass limit of the pion result given above; 
namely uj^ = 2VT 3 /(7i 2 k 4 ). An additional factor of 2 appears for the spin of 
the photon. 
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